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Abstract. The quantum spin Hall effect shares many similarities (and some 
important differences) with the quantum Hall effect for the electric charge. As with 
the quantum (electric charge) Hall effect, there exists a correspondence between bulk 
and boundary physics that allows to characterize the quantum spin Hall effect in 
diverse and complementary ways. In this paper, we derive from the network model 
that encodes the quantum spin Hall effect, the so-called Z 2 network model, a Dirac 
Hamiltonian in two dimensions. In the clean limit of this Dirac Hamiltonian, we 
show that the bulk Kane-Mele Z 2 invariant is nothing but the SU(2) Wilson loop 
constructed from the SU(2) Berry connection of the occupied Dirac-Bloch single- 
particle states. In the presence of disorder, the non-linear sigma model (NLSM) 
that is derived from this Dirac Hamiltonian describes a metal-insulator transition in 
the standard two-dimensional symplcctic universality class. In particular, we show 
that the fermion doubling prevents the presence of a topological term in the NLSM 
that would change the universality class of the ordinary two-dimensional symplcctic 
metal-insulator transition. This analytical result is fully consistent with our previous 
numerical studies of the bulk critical exponents at the metal-insulator transition 
encoded by the Z 2 network model. Finally, we improve the quality and extend the 
numerical study of boundary multifractality in the Z 2 topological insulator. We show 
that the hypothesis of two-dimensional conformal invariance at the metal-insulator 
transition is verified within the accuracy of our numerical results. 



PACS numbers: 73.20.Fz, 71.70.Ej, 73.43.-f, 05.45.Df 
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1. Introduction 

Spin-orbit coupling has long been known to be essential to account for the band structure 
of semiconductors, say, semiconductors with the zink-blende crystalline structure. 
Monographs have been dedicated to reviewing the effects of the spin-orbit coupling on 
the Bloch bands of conductors and semiconductors [![]. Electronic transport properties of 
metals and semiconductors in which impurities are coupled to the conduction electrons 
by the spin-orbit coupling, i.e., when the impurities preserve the time- reversal symmetry 
but break the spin-rotation symmetry, are also well understood since the prediction 
of weak antilocalization effects ||. Hence, the prediction of the quantum spin Hall 
effect in two-dimensional semiconductors with time-reversal symmetry but a sufficiently 
strong breaking of spin-rotation symmetry is rather remarkable in view of the maturity 
of the field dedicated to the physics of semiconductors 0, f|, ||, ||. The quantum 
spin Hall effect was observed in HgTe/(Hg,Cd)Te quantum wells two years later 0. 
Even more remarkably, this rapid progress was followed by the prediction of three- 
dimensional topological insulators || || |1(| and its experimental confirmation for Bi- 



based compounds JTT], [12], [13], [lj, [15] . 

The quantum spin Hall effect, like its relative, the quantum (electric charge) Hall 
effect, can be understood either as a property of the two-dimensional bulk or as a 
property of the one-dimensional boundary. The bulk can be characterized by certain 
integrals over the Brillouin zone of Berry connections calculated from Bloch eigenstates. 
These integrals are only allowed to take discrete values and are examples of topological 
invariants from the mathematical literature. As is well known, the topological number 



v takes integer values for the quantum (electric charge) Hall effect ||16fl . By contrast 



it takes only two distinct values {v = or 1) for time-reversal invariant, Z 2 topological 



band insulators ||, ||, [H], [lTj • Because they are quantized, they cannot change under a 



CONTENTS 



3 



small continuous deformation of the Hamiltonian, including a perturbation that breaks 
translation invariance, i.e., disorder. 

The bulk topological quantum numbers are closely connected with the existence 
of stable gapless edge states along the boundary of a topological insulator, or more 
precisely along the interface between two insulators with different topological numbers. 
The number of gapless edge modes is determined by the difference of the topological 
numbers. On the edge of a two-dimensional Z 2 topological band insulator with v — 1, 
there exists helical edge states, a Kramers' pair of counter propagating modes, which 
interpolates between the bulk valence band and the bulk conduction band. If one 
changes the Fermi energy from the center of the band gap to lower energies through 
the conduction band, one should observe a transition from a Z 2 topological insulator 
to a metal, and then from a metal to a trivial band insulator [y = 0) without helical 
edge states. Since both helical edge states and a metallic phase are stable against 
(weak) disorder (due to the quantized topological number and to weak anti-localization, 
respectively), the same sequence of phases should appear as the Fermi energy is varied 



even in the presence of disorder, as confirmed recently by numerical simulations [18, 19 



A question one can naturally ask is then whether there is any difference between the 
critical phenomena at the metal-to-Z 2 -topological-insulator transition and those at the 
metal-to-trivial-insulator transition. This is the question which we revisit in this paper, 
extending our previous studies ]19|. ^(J . It will become clear that one needs to distinguish 
between bulk and boundary properties in the universal critical phenomena. 

For the quantum (electric charge) Hall effect, the Chalker-Coddington network 
model serves as a standard model for studying critical properties at Anderson transition 
between different quantum Hall states J2T]. The elementary object in the Chalker- 
Coddington network model is chiral edge states. These edge states are plane waves 
propagating along the links of each plaquette which represents a puddle of a quantum 
Hall droplet formed in the presence of spatially slowly varying potential. They are 
chiral as they represent the mode propagating along equipotential lines in the direction 
determined by the external magnetic field. The Chalker-Coddington network model is 
a unitary scattering matrix that scales in size with the number of links defining the 
network, and with a deterministic parameter that quantifies the relative probability for 
an incoming mode to scatter into a link rotated by +ir/2 or —tt/2. By tuning this 
parameter through the value 1/2, one can go through a transition from one insulating 
phase to another insulating phase, with the topological number v changed by one. 
This remains true even when the phase of an edge state along any link is taken to 
be an independent random number to mimic the effects of static local disorder. The 
Chalker-Coddington model is a powerful tool to characterize the effects of static disorder 
on the direct transition between two successive integer quantum Hall states. It has 
demonstrated that this transition is continuous and several critical exponents at this 
transition have been measured from the Chalker-Coddington model [EI], [22] . 

The present authors have constructed in |T3 a generalization of the Chalker- 



Coddington model that describes the physics of the two-dimensional quantum spin Hall 
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effect. We shall call this network model the Z 2 network model, which will be briefly 
reviewed in section 2. As with the Chalker-Coddington model, edge states propagate 
along the links of each plaquette of the square lattice. Unlike the Chalker-Coddington 
model there are two edge states per link that form a single Kramers' doublet, which 
corresponds to helical edge states moving along a puddle of a quantum spin Hall droplet. 
Kramers' doublets undergo the most general unitary scattering compatible with time- 
reversal symmetry at the nodes of the square lattice. The Z 2 network model is thus 
a unitary scattering matrix that scales in size with the number of links defining the 
network and that preserves time-reversal symmetry. The Z 2 network model supports 
one metallic phase and two insulating phases, as we discussed earliei(||. The metallic 
phase prevents any direct transition between the insulating phases and the continuous 
phase transition between the metallic and any of the insulating phases belongs to the 
two-dimensional symplectic universality class of Anderson localization ||. 

Numerical simulations have shown that bulk properties at metal-insulator transition 
in the Z 2 network model are the same as those at conventional metal-insulator transitions 



in the two-dimensional symplectic symmetry class |19], ^0|. In fact, one can understand 
this result from the following general argument based on universality. The non-linear 
sigma model (NLSM) description is a very powerful, standard theoretical approach to 



Anderson metal-insulator transition |2"5 |. A NLSM can have a topological term if the 
homotopy group of the target manifold, which is determined by the symmetry of the 
system at hand, is nontrivial. Interestingly, in the case of the symplectic symmetry class, 
as is called the statistical ensemble of systems (including quantum spin Hall systems) 
that are invariant under time reversal but are not invariant under SU(2) spin rotation, 



the NLSM admits a Z 2 topological term |2|, [25[ [2§]. Moreover, the NLSM in the 
symplectic symmetry class with a Z 2 topological term cannot support an insulating 
phase. This can be seen from the fact that this NLSM describes surface Dirac fermions 
of a three-dimensional Z 2 topological insulator which are topologically protected from 
Anderson localization |27], |29[ . This in turn implies that any two-dimensional metal- 
insulator transition in time-reversal-invariant but spin-rotation-noninvariant systems 
should be in the same and unique universality class that is encoded by the NLSM 
without a topological term in the (ordinary) symplectic class. 

Whereas bulk critical properties at the transition between a metal and a Z 2 
topological insulator do not depend on the topological nature of the insulating phase, 
there are boundary properties that can distinguish between a topologically trivial and 
non-trivial insulating phases. Boundary multifractality is a very convenient tool to 
probe any discrepancy between universal bulk and boundary properties at Anderson 
transition [j30|, pH]]. To probe this difference, the present authors performed a multifractal 
analysis of the edge states that propagate from one end to the other in a network model 



at criticality with open boundary condition in the transverse direction [20]. It was found 



that boundary multifractal exponents are sensitive to the presence or absence of a helical 

% The presence or absence of a single helical edge state in an insulating phase is solely dependent on 
the boundary conditions which one imposes on the network model. 
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Kramers' doublet propagating along the boundary. 
The goal of this paper is 2-fold: 

(i) to establish a direct connection between the Z 2 network model and a Hamiltonian 
description of the Z 2 topological insulator perturbed by time-reversal symmetric 
local static disorder. 

(ii) to improve the quality and extend the numerical study of boundary multifractality 
in the Z 2 topological insulator. 

For item (@), in section 3, we are going to relate the Z 2 network model to a problem 
of Anderson localization in the two-dimensional symplectic universality class that is 
encoded by a stationary 4x4 Dirac Hamiltonian perturbed by static disorder that 
preserves time-reversal symmetry but breaks spin-rotation symmetry. This result is a 
natural generalization of the fact that the Chalker-Coddington network model can be 
related [32] to a 2 x 2 Dirac Hamiltonian with static disorder [p3]. In the clean limit, 
we shall characterize the Z 2 insulating phases in the 4x4 Dirac Hamiltonian by a 
Z 2 topological invariant. In particular, we show that an SU(2) Wilson loop of Berry 
connection of Bloch wave functions is equivalent to the Z 2 index introduced by Kane 
and Mele [|J. The 4x4 Dirac Hamiltonian will allow us to make contact between the 
Z 2 network model and the NLSM description of two-dimensional Anderson localization 
in the symplectic universality class derived 30 years ago by Hikami et al. in 0. In 
our opinion, this should remove any lingering doubts that the metal-insulator transition 
between a two-dimensional metallic state and a two-dimensional Z 2 insulator that is 
driven by static disorder is anything but conventional. 

For item (0), besides improving the accuracy of the critical exponents for one- 
dimensional boundary multifractality in the Z 2 network model, we compute critical 
exponents for two zero-dimensional boundaries (corners) in section 4. We shall use 
these critical exponents to verify the hypothesis that conformal invariance holds at the 
metal-insulator transition and imposes relations between lower-dimensional boundary 
critical exponents. 



2. Definition of the Z 2 network model for the quantum spin Hall effect 

The Z 2 network model is defined as follows. First, one draws a set of corner sharing 
square plaquettes on the two-dimensional Cartesian plane. Each edge of a plaquette is 
assigned two opposite directed links. This is the network. There are two types S and 
S' of shared corners, which we shall call the nodes of the network. Second, we assign to 
each directed link an amplitude ip, i.e., a complex number ip G C. Any amplitude ip is 
either an incoming or outgoing plane wave that undergoes a unitary scattering process 
at a node. We also assign a 4 x 4 unitary matrix S to each node of the network. The set 
of all directed links obeying the condition that they are either the incoming or outgoing 
plane waves of the set of all nodal unitary scattering matrices defines a solution to the 
Z 2 network model. 
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Figure 1. (a) The Z 2 network model. The solid and dashed lines represent the links 
for up and down spin electrons, respectively. The electrons are unitarily scattered at 
the nodes S and S'. The choice for the scattering basis at the nodes S and S' is shown 
in (b) and (c), respectively. 



To construct an explicit representation of the Z 2 network model, the center of each 
plaquette is assigned the coordinate (x, y) with x and y taking integer values, as is 
done in figure [I]. We then label the 8 directed links ip ncT {x,y) of any given plaquette by 
the coordinate (x,y) of the plaquette, the side n = 1,2,3,4 of the plaquette with the 
convention shown in figure [I], and the spin index o —\ or a =| if the link is directed 
counterclockwise or clockwise, respectively, relative to the center of the plaquette. The 
4x4 unitary S'-matrix is then given by 



( 



^( x ,y) 

V> 4T + 1,3/- 
\ ^(z + 1,2/- 
at any node of type S or as 

/ ^3t( x + 1 '2/- 

at any node of type S', with 

S = U(x,y)S V(x,y) 
Here, the 4x4 unitary matrix 

tQ 



\ 



S 



3?V 

^{x^y) 

n x + l,y- 



\ 



X 



(1) 



i,y-i) / 



S' 



^4T^ I - 1) y + 1) 

^( x ^y) 
^n( x ,y) 



(2) 



S' = U'(x,y)S V'(x,y). 



S :-- 



TSr 



—tQ? rs 

is presented with the help of the unit 2x2 matrix s and of the 2x2 matrix 

cos 9 sin 6 



Q := s 1 sin 9 + s 3 cos 9 



sin 9 — cos 9 



(3) 



(4) 



(5) 
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(s 1; s 2 , and s 3 are the 2x2 Pauli matrices) that are both acting on the spin indices 
<7 =T, !, together with the real- valued parameters 

r := tanhX, t := r - ^? (6) 

cosh A 

with 

{(X,0) |0 < X < oo, 0<6<7r/2}. (7) 

For later use, we shall also introduce the real-valued parameter (3 G [0, 7r] through 

r = cos/3, £ = sin/3. (8) 

The parameter 6 controls the probability of spin-flip scattering, sin 2 6. The unitary 
matrices U, V, U', V are defined as 

U(x,y) = diag(e i ^ (x ' y) ,e ix 3^), e i ^ (a;+1 ' y - 1) ,e ix i (:c+1 ' !/ - 1) ), (9a) 

V(x,y) = diag(e ix 3^'f),e ix 2 (a; ' 3/) ,e ix i (x+1 ' y - 1) ,e ix 4( :c + 1 ^- 1 )), (9b) 

U'(x,y) = diag(e ix 3^+ 1 ^+ 1 ), e ix 4( :c + 1 ^+ 1 ),e ix i (a; ' y) ,e ix 2(^)) ! (9c) 

V'{x,y) = diag(e ix ^ a;+1 ' y+1 ),e ix 3( a;+1 ' 3/+1 ),e ix 2( a; ^,e ix i (:E ' J/) ), (9d) 

where 2x n (x,y) equals a (random) phase that wave functions acquire when propagating 
along the edge n of the plaquette centered at (x, y). 

The Z 2 network model is uniquely defined from the scattering matrices S and S' . 
By construction, the S'-matrix is time- reversal symmetric, i.e., 

iS2 V ( ~ iS2 ° (10) 

is 2 J \ -is 2 J K J 

and a similar relation holds for S'. 

In [J9| , we obtained the phase diagram of the Z 2 network model shown schematically 
in figure 0(a). Thereto, (X, 6) are spatially uniform deterministic parameters that can 
be changed continuously. On the other hand, the phases \ n of all link plane waves in 
the Z 2 network model are taken to be independently and uniformly distributed random 
variables over the range [0, 27r). The line 9 = is special in that the Z 2 network model 
reduces to two decoupled Chalker-Coddington network models [19|. Along the line 
= 0, the point 

X cc = 1x1(1 + ^2)^/3 = j (11) 

realizes a quantum critical point that separates two insulating phases differing by one 
gapless edge state or, equivalently, by one unit in the Hall conductivity, per spin. 
Alternatively, 8 can also be chosen to be randomly and independently distributed at 
each node with the probability sin(2#) over the range (0, vr/2). This leaves X as the sole 
deterministic parameter that controls the phase diagram as shown in figure 0(b). When 
performing numerically a scaling analysis with the size of the Z 2 network model, one 
must account for the deviations away from one-parameter scaling induced by irrelevant 
operators. The Z 2 network model with a randomly distributed 9 minimizes such finite- 
size effects (see jl9| ). 
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(b) 

Insulator Metal Insulator 
• 1 1 ► X 

x s Xl 



" 1 
X 

Figure 2. (a) Schematic phase diagram from the analysis of the Z 2 network model 
with the constant X and 8. The metallic phase is surrounded by the two insulating 
phases with the critical points X s and Xi(> X s )) for < 8 < ir/2. The fixed point 
denoted by a filled (green) square along the boundary 8 = is the unstable quantum 
critical point located at X cc = ln(l + \/2) separating two insulating phases in the 
Chalker-Coddington model. The fixed point denoted by the filled (blue) rhombus at 
the upper left corner is the unstable metallic phase. The shape of the metallic phase is 
controlled by the symmetry crossover between the unitary and symplectic symmetry 
classes, (b) The phase diagram for Z 2 network model with randomly distributed 8 
over the range (0, w/2). 

3. Two-dimensional Dirac Hamiltonian from the Z 2 network model 

The Chalker-Coddington model is related to the two-dimensional Dirac Hamiltonian as 
was shown by Ho and Chalker in [BS]. We are going to establish the counterpart of 
this connection for the Z 2 network model. A unitary matrix is the exponential of a 
Hermitian matrix. Hence, our strategy to construct a Hamiltonian from the Z 2 network 
model is going to be to view the unitary scattering matrix of the Z 2 network model as a 
unitary time evolution whose infinitesimal generator is the seeked Hamiltonian. To this 
end, we proceed in two steps in order to present the Z 2 network model into a form in 
which it is readily interpreted as a unitary time evolution. First, we change the choice of 
the basis for the scattering states and select the proper unit of time. We then perform 
a continuum approximation, by which the Z 2 network model is linearized, so to say. 
This will yield an irreducible 4-dimensional representation of the Dirac Hamiltonian in 
(2 + l)-dimensional space and time, a signature of the fermion doubling when deriving a 
continuum Dirac Hamiltonian from a time-reversal symmetric and local two-dimensional 
lattice model. 

3.1. Change of the basis for the scattering states and one-step time evolution 

Our goal is to reformulate the Z 2 network model defined in Sec. fjin such a way that the 
scattering matrix maps incoming states into outgoing states sharing the same internal 
and space labels but a different "time" label. This involves a change of basis for the 
scattering states and an "enlargement" of the Hilbert space spanned by the scattering 
states. The parameter 6 is assumed to be spatially uniform. We choose the plaquette 
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(x, y) of the network. 

At node S of the plaquette (x, y) , we make the basis transformation and write the 
S'-matrix (Jl]) in the form 



^31 
V ^4 T J 

where we have defined 
/ 



^3T 
V^4| J 



M S =UM S U, 



(12) 



A4 



V 






-tt x _t\ 


sin# 


tt x _t y + cos 6 


r 


t t\t y _ sin 9 







r 


—t t X jt y _ cos 


t t\t y _ cos 9 


r 







tt x + t y _ sine 


r 


-tt x _t y 


cos# 


-t t x j% sine 






y 



and 



(13) 



(14) 



U{x,y) = diag(e ix i (x ' y) ,e ix 3( a; ^,e ix 2( a; ^,e i ^^^). 
Here given n = 1, 2, 3, 4 and a =f, I, we have introduced the shift operators acting on 

t±1Pn*( X >y) '■= ^na{%± 1,2/), ^±^n CT (^,2/) : = ^n(^,y ± !), (15) 

and similarly on the phases Xni x iU) ^ [0,2tt). We note that the scattering matrix A/s 
is multiplied by the unitary matrix IA from the left and the right in (0), because the 
Kramers' doublet acquires exactly the same phase Xn when traversing on the edge n of 
the plaquette (x, y) before and after experiencing the scattering As at the node S. 

At node S' of the plaquette (x,y), we make the basis transformation and rewrite 
the scattering matrix S' (|2|) into the form 

M s , = UN S ,U, 



^3T 
\ J 

where we have defined 



= : M v 



V ^4 T ) 



(16) 



A/;, 



/ 





—t sin e r 


-tt x + t\ 


cose \ 




t t x _t y _ sin 9 


o t t x _ti cos e 


r 






r 


tt x + t\ cose o 


-tt\t\ 


sin e 


V 


-t t x _t y _ cose 


r tt x _t y _sine 





) 



;i7) 



As it should be 

M\M S = M j sl M s , = 1. (18) 

Next, we introduce the discrete time variable / G Z as follows. We define the 
elementary discrete unitary time evolution to be 

M s 
M s , 



v^-i y 



v^-i y 



(19) 
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Here, to treat on equal footing the nodes of type S and S', we have enlarged the scattering 
basis with the introduction of the doublets 



0H 



(20) 



Due to the off-diagonal block structure in the elementary time evolution, it is more 
convenient to consider the "one-step" time evolution operator defined by 



{4+1 \ 



M S M S , 




1+2 



{^+1 \ 

\4-l J 



M 



SS' 



{ 4 +l \ 
\4>-i J 



The two Hamiltonians generating this unitary time evolution are then 



n 



SS' 



-ilnM 



SS'' 



^S'S 



+ilnA<c, 



S'S- 



(21) 



(22) 
(23) 



Evidently, the additivity of the logarithm of a product implies that 
^SS' = ^S's- 

From now on, we will consider H ss , exclusively since M s , s = exp(i7i s , s ) merely 
duplicates the information contains in M ss , = exp(i7-^ ss ,). 



3.2. Dirac Hamiltonian close to 6 = 

In this section, we are going to extract from the unitary time-evolution (l21~l)-(|23"D of 
the Z 2 network model a 4 x 4 continuum Dirac Hamiltonian in the close vicinity of the 
quantum critical point 

(M) CC :=(0,7r/4). (24) 

To this end and following ||32|| , it is convenient to measure the link phases \ n in = 
1,2,3,4) relative to their values when they carry a flux of tt per plaquette. Hence, we 
redefine 



71 

X4+2 



on all plaquettes. 

Our strategy consists in performing an expansion of 

^ss' = +i ln -^ss' = +i ( m -Ms + In -M s ) = 
defined in ( p2|) to leading order in powers of 



ilnX S / S = H v 



S'S 



+ ' 



(25) 

(26) 
(27) 
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with n = 1, 2, 3, 4 where d x>y is the generator of infinitesimal translation on the network 
(the two-dimensional momentum operator). 

When 9 = 0, the unitary time-evolution operator at the plaquette (x, y) is given by 



^+1 




(2? 



1+2 

(o) _ / AW» 

^ss' - £(0) C (0) I > 



whereby 



M®=M§W°\ (30) 
M®-( ° Ai0) ) M^-( ° (31) 
with the 2x2 operator-valued matrices 

1 e i(x 3 +x 2 ) cos p -ie^Uple 1 ^ sin (3 J ' 

(0) _ /e ix 2^e ix isin/3 e^+^cos/? \ 
1 ]jbu+Xi)coaP -ie ix H"Lt y + e ix * sin f3 )' 

(0) _ / e^i+XaJcos/S -ie ix H x + t + e ix * sin f3 \ 

\e ix H x _t y _e ix zsin(3 ie^+xj CO s /3 )' ( t 

d(0) _ / e^+Xi) C os/3 e^^e^sin/A 
I -ie^rtV** sin/3 ie^+Xa) CO s/3 J" 

Observe that in the limit = 0, the Z 2 network model reduces to two decoupled U(l) 
network models where each time evolution is essentially the same as the one for the 



U(l) network model derived in p2 



In the vicinity of the Chalker-Coddington quantum critical point ([24]), we find the 
4x4 block diagonal Hamiltonian 

Hi = ( D + I ) (36) 

where the 2x2 block are expressed in terms of linear combinations of the 2x2 unit 
matrix cx and of the Pauli matrices a x , a y , and a z according to 

D + = a z {-id x + A x ) - a x (-id y + A y ) - o y m + a A , (37) 

and 

D_ = -a y (-id x - A x ) + a z (-id y - A y ) + a x m + a A . (38) 
Thus, each 2x2 block Hamiltonian is of the Dirac form whereby the linear combinations 

A) : = ~(Xl +X2 + X3 + X4), (ActAy) : = (-Xi +X3,X2 -Xj> (39) 

enter as a scalar gauge potential and a vector gauge potential would do, respectively. 
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Any deviation of 9 from 9 = lifts the reducibility of (0). To leading order in 9 
and close to the Chalker-Coddington quantum critical point (|24|), 



M ss , = (Mf + 9Mf + ■■■) (MP + 6M$> 



e(M?M§ ) + M§ ) M$)+--- (40) 





' o 


i 









1 




— i 




(: 







with 

where we have set m 

^SS'= | ~t ^ I. D e :=9[ ; ; ) (43) 

to this order. 

Next, we perform a sequence of unitary transformation generated by 




e 



i7RT y /4 Q \ / e ~ma x /4 Q 




yielding 



aa tl_ 



with a = \f 7 l9 and 

W± = o x (-id x ± A,) + <7„ ± A y ) ± a z m + or A . (46) 

The 2x2 matrices 7i + and 7i_ describe a Dirac fermion with mass ±m in the presence of 
random vector potential ±(A„, A ) and random scalar potential A , each of which is an 
effective Hamiltonian for the plateau transition of integer quantum Hall effect [KI3, . 
The TC± sectors are coupled by the matrix element aa . 

The 4x4 continuum Dirac Hamiltonian TC can be written in the form 

W = (-id x a x - ifyo-J ® r + (A x a x + + ma z ) ® 

+ A a <g> r + a a <g> r x , (47) 

where r is a unit 2x2 matrix and r x , r y , and r z are three Pauli matrices. The 
Hamiltonian (|4"7 ) is invariant for each realization of disorder under the operation 

TH* T _1 = H, T := ia y ® r x , (48) 

that implements time-reversal for a spin-1/2 particle. 

The Dirac Hamiltonian fl4"7|) is the main result of this subsection. It is an effective 
model for the Anderson localization of quantum spin Hall systems, which belongs to 
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the symplectic class in view of the symmetry property fl48|). The Anderson transition 
in the Dirac Hamiltonian ( |4"?| ) should possess the same universal critical properties as 
those found in our numerical simulations of the Z 2 network model. In the presence of 
the "Rashba" coupling a, there should appear a metallic phase near m = which is 
surrounded by two insulating phases. In the limit a — > 0, the metallic phase should 
shrink into a critical point of the integer quantum Hall plateau transition. 

The 4x4 continuum Dirac Hamiltonian 7i should be contrasted with a 2 x 2 
Hamiltonian of a Dirac particle in random scalar potential, 

H 2 = -id x a x - id y a y + V(x, y)a , (49) 

which has the minimal dimensionality of the Clifford algebra in (2 + l)-dimensional 
space time and is invariant under time-reversal operation, cr T-Qcr = 7Y 2 . The 2x2 Dirac 
Hamiltonian ( |4"9| ) is an effective Hamiltonian for massless Dirac fermions on the surface of 
a three-dimensional Z 2 topological insulator. After averaging over the disorder potential 
V, the problem of Anderson localization of the surface Dirac fermions is reduced to 



a NLSM with a Z 2 topological term [25. Era. Interestingly, this Z 2 topological term 



prevents the surface Dirac fermions from localizing [E7I, EBJ. It is this absence of two- 



dimensional localization that defines a three-dimensional Z 2 topological insulator [ 29[| . 
In contrast, the doubling of the size of the Hamiltonian (|47|) implies that the NLSM 
describing the Anderson localization in the 4x4 Hamiltonian ([45"D does not come with 
a Z 2 topological term, because two Z 2 topological terms cancel each other. We can thus 
conclude that the critical properties of metal-insulator transitions in the Z 2 network 
model are the same as those in the standard symplectic class, in agreement with results 



of our numerical simulations of the Z 2 network model p9| , p0 



Before closing this subsection, we briefly discuss the Dirac Hamiltonian ([47|) in the 
clean limit where A = A x = A y = 0. Since the system in the absence of disorder is 
translationally invariant, momentum is a good quantum number. We thus consider the 
Hamiltonian in momentum space 

n{k)= ( K^ + k y a y + ma z «*o ) (5Q) 

where the wave number k = (k x ,k y ). When a = 0, the Hamiltonian (^) becomes a 
direct sum of 2 x 2 Dirac Hamiltonian with mass of opposite signs. This is essentially 
the same low-energy Hamiltonian as the one appearing in the quantum spin Hall effect 
in HgTe/(Hg,Cd)Te quantum wells 0. 

3.3. Z 2 topological number 

We now discuss the topological property of the time-reversal invariant insulator which 
is obtained from the effective Hamiltonian (^) of the Z 2 network model in the absence 
of disorder. The topological attribute of the band insulator is intimately tied to the 
invariance 

e -1 ft(-fc)e = H(k) (51) 
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under the operation of time-reversal represented by 



(52) 



where /C implements complex conjugation. We are going to show that this topological 
attribute takes values in Z 2 , i.e., the Z 2 index introduced by Kane and Mele Q. 

We begin with general considerations on a translation-invariant single-particle 
fermionic Hamiltonian which has single-particle eigenstates labeled by the wave vector 
k taking values in a compact manifold. This compact manifold can be the first Brillouin 
zone with the topology of a torus if the Hamiltonian is defined on a lattice and 
periodic boundary conditions are imposed, or it can be the stereographic projection 
between the momentum plane M 2 and the surface of a three-dimensional sphere if the 
Hamiltonian is defined in the continuum. We assume that (i) the antiunitary operation 
= — _1 = — 0t that implements time-reversal leaves the Hamiltonian invariant, 
(ii) there exists a spectral gap at the Fermi energy, and (iii) there are two distinct 
occupied bands with the single-particle orthonormal eigenstates \u & (k)) and energies 
E & (k) labeled by the index a — 1, 2 below the Fermi energy. All three assumptions are 
met by the 4x4 Dirac Hamiltonian (f|7]), provided that the mass m is nonvanishing. 

Because of assumptions (i) and (ii) the 2x2 unitary sewing matrix with the matrix 
elements w^k) defined by 



i.e., the overlaps between the occupied single-particle energy eigenstates with momentum 
— k and the time reversed images to the occupied single-particle energy eigenstates with 
momentum k, plays an important role [[[7]]. The matrix elements (|53|) obey 




a,b = 1 



2, (53) 



w^(k) = {u a (-k)\(e\ Ui (k))) 
= <«t(fc)|(et|u 4 (-fc)>) 
= -<« g (fc)|(6|« a (-fc))) 



= -w f J-k), a, b =1,2. (54) 



We used the fact that is antilinear to reach the second equality and that it is 
antiunitary with 2 = —1 to reach the third equality. Hence, the 2x2 unitary sewing 
matrix w(k) with the matrix elements (|53| ) can be parametrized as 




(55) 



with the three complex-valued functions 
w n (k) = -w n (-k), ? 
We observe that w(k) reduces to 



w 22 (k) 



w 22 (-k) 



w 12 (k). 



(56) 




(57) 
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A a (k) :=<u a (fc)|d«j(fc)> = (« a (fc) 



,(fc) )d* =^(fc)dfc , (58) 



for some real-valued f(k) at any time-reversal invariant wave vector k ~ —k (time- 
reversal invariant wave vectors are half a reciprocal vector for a lattice model, and or 
oo for a model in the continuum). 

As we shall shortly see, the sewing matrix imposes constraints on the U(2) 
Berry connection 

_d_ 

dic^"' / ~ ^ v 

where the summation convention over the repeated index /i is understood (we do 
not make distinction between superscript and subscript). Here, at every point k in 
momentum space, we have introduced the U(2) antihermitian gauge field A^{k) with 
the space index p, — 1, 2 and the matrix elements 

A u k )=-{ A my (59) 

labeled with the U(2) internal indices a, b = 1, 2, by performing an infinitesimal 
parametric change in the Hamiltonian. We decompose the U(2) gauge field (|58"D into 
the U(l) and the SU(2) contributions 

A^k) = al{k) P ^ + % (k).^ (60) 

where p is a 2 x 2 unit matrix and p is a 3 vector made of the Pauli matrices p x , p , 
and p z . Accordingly, 

A^ 2 \k) = A u ^(k)+A slJi2 \k). (61) 

Combining the identity 2 = — 1 with the (partial) resolution of the identity 
J2a=i2 l n a(k))( n a(k)l for the occupied energy eigenstates with momentum k yields 

£e|« a (fc))<« a (fc)|e = -i, (62) 

a=X,2 

where the proper restriction to the occupied energy eigenstates is understood for the unit 
operator on the right-hand side. Using this identity, we deduce the gauge transformation 

w u ^- k) )) 

°> / / a,S=i,2 

= -w(k)Al(k)w\k) - w{k)d^{k) 

= +w(k)Al(k)w\k) - w^d^ik) (63) 

that relates the U(2) connections at ±fe. For the U(l) and SU(2) parts of the connection, 

aJ(-fc) = aJ(fc)-25 M C(fc), (64) 
o M (-fc) ■ p = a M (fe) ■ w(k)p T w j (k) - 2i w(k)d^w j (k), (65) 

where we have decomposed w{k) into the U(l) (e 1 ^) and SU(2) (w) parts according to 

w (k) = eW*w(k), (66) 

(note that this decomposition has a global sign ambiguity, which, however, will not 
affect the following discussions). 
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Equipped with these gauge fields, we introduce the U(2) Wilson loop 



U(1)M |,,/ SU(2) 

where the U(l) Wilson loop is given by 



c 

W m] [C]xW sm] [C], (67) 



W m [C] :=exp ^ u «(fc)j , 
)p is given by 
^trPexp I^A sv{2 \k) J 



(6* 



while the SU(2) Wilson loop is given by 

^su( 2 )[C] := JtrPexp ( </ ^ u ^(fc) | . (69) 



The symbol V in the definition of the U(2) Wilson loop represents path ordering, while 
C is any closed loop in the compact momentum space. 

By construction, the U(2) Wilson loop ( |6"7] ) is invariant under the transformation 

A"(fc) -> f/ f (fc) A"(fc)l/(fc) + U\k)d"U(k) (70) 

induced by the local (in momentum space) U(2) transformation 

|u a (fc)> - \u h (k))U-Jk) (71) 

on the single-particle energy eigenstates. Similarly, the SU(2) and U(l) Wilson loops 
are invariant under any local SU(2) and U(l) gauge transformation of the Bloch wave 
functions, respectively. 

When C is invariant as a set under 

k -> -k, (72) 
the SU(2) Wilson loop W 7 ^^ [C] is quantized to the two values 

W sm [C]=±l (73) 
because of time-reversal symmetry. Furthermore, the identity 

K K 

W SV{2) [C}= n Pf(tS(JO), (74) 
Kec 

which we will prove below, follows. Here, the symbol Pf denotes the Pfaflian of an 
antisymmetric matrix, and only the subset of momenta K G C that are unchanged 
under K — > —K contribute to the SU(2) Wilson loop. According to (|54|), the sewing 
matrix at a time-reversal symmetric wave vector is an antisymmetric 2x2 matrix. 
Consequently, the SU(2) part of the sewing matrix at a time-reversal symmetric wave 
vector is a real- valued antisymmetric 2x2 matrix (i.e., it is proportional to ip y up to a 
sign). Hence, its Pfaflian is a well-defined and nonvanishing real- valued number. 

Before undertaking the proof of (0), more insights on this identity can be obtained 
if we specialize to the case when the Hamiltonian is invariant under any U(l) subgroup 
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of SU(2), e.g., the ^-component of spin a z . In this case we can choose the basis states 
which diagonalize a z \ a z \u x (k)) = +\u 1 (k)), a z \u 2 (k)) = —\u 2 (k)). Since the time- 
reversal operation changes the sign of a z , the sewing matrix takes the form 

n p-ix(fc) 
-e-«-*> 



which, in combination with flB"i| ) and (^) implies the transformation laws 

= + d , lx(k) + X (-k)} , (76) 

K(-k) = -o-(fc) + |x(fc) - x(-fe)] • (77) 

We conclude that when both the z component of the electron spin and the electron 
number are conserved, we can set 

«Z(*0 = <(k) = 0, A^\k) = aj(fc)| + a*(fc)|, (78) 

and use the transformation law 

= ^ K(-fc) + al(-k)] = A^£{k) - id vX {k). (79) 

With conservation of the z component of the electron spin in addition to that of 
the electron charge, the SU(2) Wilson loop becomes 



Wsu ( 2)[C] 



UiVe^ij A sv{2 \k)\ (80) 



itr exp I a*(fc)^ 



= oosfiy'a*(fc)dF 

We have used the fact that is traceless to reach the last line. This line integral can 
be written as the surface integral 

ja^k)^ = J d?k£^d„a z v {k) (82) 



v 



by Stokes' theorem. Here, T> is the region defined by dT> = C, and covers a half of the 
total Brillouin zone (BZ) because of the condition ([721). In turn, this surface integral is 
equal to the Chern number for up-spin fermions, 

Ch T := / ^^ D S(fc) (83) 

JBZ ^ 
d 2 k 



bz 27ri 
d 2 k 



^\k) 



v 2ti'\ 
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d 2 k 



0„ 



I'D 

— i 



2ni~ 



U(2). 



2m 



where we have used the transformation law (ITS) to deduce that 



-k) 



,U(2) 



22 



(k) 



(84) 



55) 



to reach the fourth equality. 

To summarize, when the z component of the spin is conserved, the quantized SU(2) 
Wilson loop can then be written as the parity of the spin Chern number (the Chern 
number for up-spin fermions, which is equal to minus the Chern number for down-spin 
fermions) |, |, g, 



Wa 



\Ch T 



'SU(2)[CJ — ( 

Next, we apply the master formula ([74]) to the 4x4 Dirac Hamiltonian (|5 
this end, we first replace the mass m by the fc-dependent mass 

m - Ck 2 , C>0, 



m. 



(86) 
To 

(87) 



and parametrize the wave number k as 



-oo < k < oo, 



< If < IT. 



OO. 



Without loss of generality, we may assume a > 0. The mass m k is introduced so that 
the SU(2) part of the sewing matrix is single- valued in the limit \k\ 
We then perform another series of unitary transformation with 

3 i7T(7 z /4 



u 



(Tn 







i«7. 

to rewrite the Hamiltonian 










-incr z /A 



(89) 



in the form 



H{k) := U ] K{k)U 



K°x + k y a y + ( a + im k) a o 



K a x + k y a y + ( a ~ im k) ff 



(90) 

The four eigenvalues of the Hamiltonian (|90|) are given by E(k) = ±A^, ±A^, where 

Af = y/(k±a)* + ml. (91) 
The occupied eigenstate with the energy Ei(k) = — A^ reads 



l -K \ 



|iti(<p,fc)) 



2X1 



-k + a + im k 



a — imje 



(92) 
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Table 1. 6f at the time -reversal invariant momenta k = and k = ±oo, when 
m — Ca 2 < (a) and m — Co? > (b). It is assumed that < arctan(|m|/a) < ir/2. 



a) m — 


Ca 2 


< 






u 

h 




— oo 


U 


+oo 


a+ 

k 




/o 
7T/2 


arctan(— m/a) 


/o 
7T/2 


"k 




tt Id 


djICLdlll — 1 11 JUc) — /I 


TT /9 






i 


i 


i 


(b) m - 


-Ca 


2 > 






k 




— oo 





+oo 






-3tt/2 


— arctan(m/a) 


tt/2 


K 




3tt/2 


7r — arctan(m/o;) 


-tt/2 






i 


— i 


i 



and the occupied eigenstate with the energy E 2 (k) = — A^ is 

|u 2 (<M)) 



2A; 



-A+e"^ 
k + a + im. 



\ (k + a + iffi/Je" 
Notice that (^(y?, &)) = |wi(y + — 



ly 



(93) 



The 2x2 sewing matrix is obtained from the eigenstates (|92|)-(p3|) as 

w(<p,k) := ((u & (<p,-k)\Q\ui(<p,k)) 



a,6=i,2 



/ 



, — (A; - ct + imj 
V Ar v 



which is decomposed into the U(l) part, 
iexp (i^ + i(0+ + 0*)/2), 
and the SU(2) part, 



(94) 



(95) 



(96) 



(97) 



For the SU(2) sewing matrix (|9"6|), there are two momenta which are invariant under 
inversion k — > —k, namely the south K = and north K = oo poles of the stereographic 
sphere. The values of &t at these time-reversal momenta are listed in table [[[ The 



w(k) 



j e i(^ + -^)/2 

i e -i«-0fc )/2 o 



of the sewing matrix. Here, we have defined 9^ through the relation 
e*? = -L[k±{a-im k )]. 
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QSH 

(b) 




Figure 3. (a) Quantum spin Hall droplet immersed in the reference vacuum [in real 
space (x,y) € R 2 ]. (b) The Z 2 network model or its tight-binding equivalent when 
x < is separated from the reference vacuum at x > by the vertical boundary x = 
[in real space (x, y) € R 2 ]. 



Pfaffian of the sewing matrix at the south and north poles of the stereographic sphere 
are 

Pf w(0) = - sg n(m - Ca 2 )Pf (-ip y ), (98) 

PftZ;(oo)=Pf(-ip y ), (99) 

respectively. Hence, 

W sv{2) [C] = - S gn(m) (100) 

for any time-reversal invariant path C passing through the south and north poles, where 
we have suppressed Ca 2 by taking the limit Ca 2 /\m\ — > 0. 

The value ( |100|) taken by the SU(2) Wilson loop thus appears to be ambiguous since 
it depends on the sign of the mass m. This ambiguity is a mere reflection of the fact that, 
as noted in 0], the topological nature of the Z 2 network model is itself defined relative 
to that of some reference vacuum. Indeed, for any given choice of the parameters (X, 9) 
from figure that defines uniquely the bulk properties of the insulating phase in the Z 2 
network model, the choice of boundary conditions determines if a single helical Kramers' 
doublet edge state is or is not present at the boundary of the Z 2 network model. In view 
of this, it is useful to reinterpret the Z 2 network model with a boundary as realizing a 
quantum spin Hall droplet immersed in a reference vacuum as is depicted in figure 0(a). 
If so, choosing the boundary condition is equivalent to fixing the topological attribute of 
the reference vacuum relative to that of the Z 2 network model, for the reference vacuum 
in which the quantum spin Hall droplet is immersed also has either a trivial or non- 
trivial Z 2 quantum topology. A single helical Kramers' doublet propagating unhindered 
along the boundary between the quantum spin Hall droplet and the reference vacuum 
appears if and only if the Z 2 topological quantum numbers in the droplet and in the 
reference vacuum differ. 

In the low-energy continuum limit (|50"D, a boundary in real space can be introduced 
by breaking translation invariance along the vertical line x = in the real space 
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Figure 4. Momentum space k y ) £ M. 2 is discretized with the help of a rectangular 
grid on which two paths are depicted. The red path that is restricted to the upper left 
quadrant is not invariant as a set under the inversion (k x , k y ) — * —(k x , k y ). The blue 
path with its center of mass at the origin is. This path is assembled out of 16 links: 

(*0>H) = — (*15'*o)> (h>H) = — (*14>*15); fe^'s) = — (*13>*14)i (*3) = — (h2)*13)i 

(241*5) = — (*ii)*i2)j (*5>*e) = — (*io>*ii)) (*6)*V) = — (*9>*io)j (*7) *s) = — (*8i*9)- Sites 
i = i 8 = i 16 along the path are the only ones invariant under (fc^., fc ) — ► — (k x , k y ). 



(x, y) G IR 2 through the profile [see figure H(b)] 

— m, if x — > — oo, 

m(x,y) — m(x) = { (101) 

+m, if x — > +oo, 

for the mass. 

We close Sec. ||| with a justification of the master formula ([T4]). To this end, we 
regularize the continuum gauge theory by discretizing momentum space (figure (|). We 
use the momentum coordinate i 6 Z 2 on a rectangular grid with the two lattice spacings 
Ak^ > 0. To each link from the site i to the nearest-neighbor site i + fi of the grid, we 
assign the SU(2) unitary matrix 

^seV, (102) 

which is obtained by discarding U(l) part of the U(2) Berry connection. Consistency 
demands that 

Ui+fri — U iti+/i < — > A i+ ^ = A^ i+ ^. (103) 
We define the SU(2) Wilson loop to be 

W sv{2) (i , - , == \^{U h)H U HH ■ ■ ■ U iN i , Q ) (104) 
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where i n and i n+l are nearest neighbors, i.e., their difference i n+1 — i n = r\ n is a unit 
vector rj n . The Wilson loop is invariant under any local gauge transformation by which 

tW - ViUi,i+»V i+li (105) 

where the V^'s are U(2) matrices. Observe that the cyclicity of the trace allows us to 
write 

W sv(2) (i ,---,i N -i) = otrfc/, , , r ■■■U i ^ A U inA U ii ---U i (106) 



'JV-1''0 '0>'l '1>'2 
"2" "T^ "T - 1 - ~2~ 



To make contact with the master formula (0), we assume that the closed path with 
vertices ip parametrized by the index £ = 0, N — 1 obeys the condition that 



N—n 



is the wave vector — Y^n=i Vm^k^ 1 " 



%-i i s the wave vector —r] 1 Ak fJ ' 1 , 

1 is the wave vector 0, (107) 

1 1 is the wave vector +rj 1 Ak^ 1 , 

i n is the wave vector + J2m=i Vm^k^ , 

with r\ m = ±1 for m = l,---,N/2 in order to mimic after discretization the 
condition that the closed path entering the Wilson loop is invariant as a set under 
the inversion (|72|) . 

On the discretized momentum lattice the sewing matrix ( p3[) is defined by 

Mtf :=<«aH)l©l«6(0>> ( 108 ) 
which obeys the condition 

w_ { = -wj, (109) 

i.e., the counterpart to the relation (|54|). This implies that w io and u; i are 
antisymmetric unitary 2x2 matrices. Furthermore, the sewing matrix w i must also 
obey the counterpart to (|63|), namely 

i ■ ,. , = «•/'.>;. (HO) 

It now follows from (|107|) and ( |110|) that 



^ i =Wi Uj i w\ , 



^ i =w t U?i w\ , (111) 

l N-l-n' l N-n 'n+l l n' l n + l l n v ' 



ft 



T T+ 1 



N 
T 



N 
T " 



T ~T " 
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In particular, we observe that 

t U i i =w { Uji w\ U { i =w i w\ U} t Ui i =w i w\ , (112) 

since is a 2 x 2 antisymmetric unitary matrix, i.e., is the second Pauli matrix 
up to a phase factor, while 

(p- n) T p 2 = -p 2 (p- n) (113) 

holds for any three- vector n contracted with the three- vector p made of the three Pauli 
matrices. By repeating the same exercise a second time, 

l N-2' l N-l \ S JV-1''0 'Oi'l/ '1>'2 '2 s l>'2 l l \ H '0 / 1>2 

one convinces oneself that the dependences on the gauge fields A ,• and A , A,- ,• 
and A,- ,■ , and so on until A- and A ,■ at the level n of this iteration 

t JV-2''iV-l 'n-1 ' 'n t iy -n' t JV-n+l 

cancel pairwise due to the conditions (|107|) — ( JTTG[ ) implementing time-reversal invariance. 
This iteration stops when n = N/2, in which case the SU(2) Wilson loop is indeed solely 
controlled by the sewing matrix at the time-reversal invariant momenta corresponding 



to t = and £ = N/2, 

w sv(2)(h> • • • > = l tT ( w i N/2 w i) ■ ( 115 ) 

Since i and i N , 2 are invariant under momentum inversion or, equivalently, time- 
reversal invariant, 

w lN/2 =e ia ^ip 2 , w lQ =e ia °ip 2 (116) 
with a N / 2 i °-q — 0, 7r. Here, the Z 2 phases e 1Qjv / 2 and e ia ° are none other than the Pfaffians 



e B */» = PfK , e-o = P f^ o ), (H7) 



respectively. Hence, 

W SU(2)(*0> •••'%-!) = Pi { W i N/2 ) ^2 X Pf ( W ! ) (-^2) 

= pf KJ pf K) (118) 

is a special case of (|74]) . (Recall that 1/^ and are real- valued.) 



4. Numerical study of boundary multifractality in the Z 2 network model 

In [20], we have shown that (i) multifractal scaling holds near the boundary of the 



Z 2 network model at the transition between the metallic phase and the Z 2 topological 
insulating phase shown in figure 0, (ii) it is different from that in the ordinary symplectic 
class, while (iii) bulk properties, such as the critical exponents for the divergence of 
the localization length and multifractal scaling in the bulk, are the same as those in 
the conventional two-dimensional symplectic universality class of Anderson localization. 
This implies that the boundary critical properties are affected by the presence of the 
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helical edge states in the topological insulating phase adjacent to the critical point. 
In this work, we improve the precision for the estimate of the boundary multifractal 
critical exponents. We also compute numerically additional critical exponents that 
encode corner (zero-dimensional) multifractality at the metal-to-Z 2 -topological-insulator 
transition. We thereby support the claim that conformal invariance is present at 
the metal-to-Z 2 -topological-insulator transition by verifying that conformal relations 
between critical exponents at these boundaries hold. 



4-1. Boundary and corner multifractality 

To characterize multifractal scaling at the metal-insulator transition in the Z 2 network 
model, we start from the time-evolution of the plane waves along the links of the network 
with the scattering matrices defined in (0)-(0) at the nodes S and S'. To minimize finite 
size effects, the parameter 9 in (|T])-(|3) is chosen to be a random variable as explained 
in Sec. [2} We focus on the metal-insulator transition at X = X t = 0.971 as shown in 
figure 0(b). 

When we impose reflecting boundary conditions, a node on the boundary reduces 
to a unit 2x2 matrix. When the horizontal reflecting boundaries are located at nodes 
of type S', as shown in figure |(a), there exists a single helical edge states for X > X l . 
The insulating phase X > X { is thus topologically nontrivial. 

For each realization of the disorder, we numerically diagonalize the one-step time- 
evolution operator of the Z 2 network model and retain the normalized wave function 
ip a (x, y), after coarse graining over the 4 edges of the plaquette located at (2, y), whose 
eigenvalue is the closest to 1. The wave function at criticality is observed to display the 
power-law dependence on the linear dimension L of the system, 

£ \^(x,y)\ 2 ^L~ A ^-^. (119) 

<r=U 

The anomalous dimension Ag° , if it displays a nonlinear dependence on g, is the 
signature of multifractal scaling. The index ( indicates whether the multifractal scaling 
applies to the bulk (( = 2), the one-dimensional boundary (( = 1), or to the zero- 
dimensional boundary (corner) (( = 0), provided the plaquette (x, y) is restricted 
to the corresponding regions of the Z 2 network model. For ( = 1 and 0, the 
index v distinguishes the case v = O when the ^-dimensional boundary has no edge 
states in the insulating phase adjacent to the critical point, from the case v = Z 2 
when the ^-dimensional boundary has helical edge states in the adjacent insulating 
phase. We ignore this distinction for multifractal scaling of the bulk wave functions, 

A (2,0) = A (2,Z 2 ) = A ( 2)j 

since bulk properties are insensitive to boundary effects. We 
will also consider the case of mixed boundary condition for which we reserve the notation 
v = Z 2 |0. 

It was shown in [|3T| that boundary multifractality is related to corner multifractality 
if it is assumed that conformal invariance holds at the metal-insulator transition in the 
two-dimensional symplectic universality class. Conversely, the numerical verification of 
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Figure 5. (a) Boundary multifractality is calculated from the wave function 
amplitudes near a one-dimensional boundary. Periodic (reflecting) boundary 
conditions are imposed for the horizontal (vertical) boundaries. (b) Corner 
multifractality is calculated from the wave function amplitudes near a corner with 
the wedge angle •& = it/2. Reflecting boundary conditions are imposed along both 
vertical and horizontal directions. The relationship between the scattering matrix at 
a node of type S' and the scattering matrix at a node of type S implies that it is a 
vertical boundary located at nodes of type S that induces an helical edge state when 
X >X,. 



this relationship between boundary and corner multifractality supports the claim that 
the critical scaling behavior at this metal-insulator transition is conformal. So we want 
to verify numerically if the consequence of the conformal map w = z®^ , namely 

Af^ = ^A^ (120) 

where d is the wedge angle at the corner, holds. Equivalently, f^' u \ce), which is defined 
to be the Legendre transformation of A^g' u ^ + dq, i.e., 

/«'"»(a„) = - A<<» - dq + C, (122) 



must obey 



^-d = -M' v) -d), (123) 



a 



/ (M («) = ^[/ (lll/) («)-l], (124) 

if conformal invariance is a property of the metal-insulator transition. 

To verify numerically the formulas ( |120| ), ( |123| ), and ( |124| ), we consider the Z 2 
network model with the geometries shown in figure ||. We have calculated wave functions 
for systems with the linear sizes L = 50, 80, 120, 150, and 180 for the two geometries 
displayed in figure ||. Here, L counts the number of nodes of the same type along a 
boundary. The number of realizations of the static disorder is 10 5 for each system size. 
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Figure 6. (a) The boundary (filled circles, red) and corner with 8 = it/2 (open circles, 
blue) anomalous dimensions at the mctal-to-Z 2 -topological-insulator transition. The 
solid curve is computed from (120) by using the boundary anomalous dimension as an 
input. The rescalcd confirming the reciprocal relation for boundary and corner 

multifractality arc shown by upper (magenta) and lower (green) triangles, respectively, 
(b) The multifractal spectra for the boundary (filled circles, red) and the corner (open 
circles, blue). The solid curve is computed from (123) and (124). 



Figure |6](a) shows the boundary anomalous dimensions Aq (filled circles) and 
the corner anomalous dimensions A g ' Z ^ (open circles). In addition, the anomalous 
dimensions A^^ are shown by upper and lower triangles for boundary and corner 
anomalous dimensions, respectively. They fulfill the reciprocal relation 

Af "> = (125) 

derived analytically in ||34|| . Since the triangles and circles are consistent within error 



bars, our numerical results are reliable, especially between < q < 1. If we use 
the numerical values of A q 1 ' Z ^ as inputs in ( [120| ) with = ir/2, there follows the 
corner multifractal scaling exponents that are plotted by the solid curve. Since the 
curve overlaps with the direct numerical computation of Aq ' 2 ^ within the error bars, 
we conclude that the relation (|120|) is valid at the metal-to-Z 2 -topological-insulator 
transition. 

Figure ||(b) shows the boundary (filled circles) and corner (open circles) multifractal 
spectra. These multifractal spectra are calculated by using (|119|) , ( |121| ), and (|122|) . The 
numerical values of a§ are 

a (i,z 2 ) = 2 Q91 ± o 0Q2 ; (126) 

ajO'Zy = 2.179 ±0.01. 

The value of a^'^ is consistent with that reported in p0[ , while its accuracy is 
improved. The solid curve obtained from the relations ( |123|) and ( |124|) by using p 1 ^ (a) 
as an input, coincides with f^ ' Z2 \a). We conclude that the hypothesis of conformal 
invariance at the quantum critical point of metal-to-Z 2 -topological-insulator transition 
is consistent with our numerical study of multifractal scaling. 
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Figure 7. (a) The z dependence of (m|\I r | 2 ) 2L at the mctal-to-Z 2 -topological- 
insulator transition in the cylindrical geometry for L = 50, 80, 120, 150, 180 from the 
top to the bottom, (b) The z dependence of a (z) (»)and c(z) (■) extrapolated 
from the system size dependence of (ln|\&| 2 ) z L averaged over a small interval of z's. 
ao(z) and c(z) at z = 0, 1 without averaging over a small interval of z's are shown 
by open circles and squares, respectively. The solid line represents the bulk value of 
off = 2.173 computed in [[si). The asymmetry with respect to z = 0.5 is due to 
statistical fluctuations. 



At last, we would like to comment on the dependence on z of 



M*IV = ^5>( E \^,y)\ 2 ) (128) 

found in [20]. Here, z = (x — 1)/2L, while x and y denote the positions on the 
network along its axis and along its circumference, respectively (our choice of periodic 
boundary conditions imposes a cylindrical geometry). The overline denotes averaging 
over disorder. Figure |7|(a) shows the z dependence of (In |\I/| 2 ) Z L for different values of 
L in this cylindrical geometry at the metal-to-Z 2 -topological-insulator transition. We 
observe that (In |^/| 2 ) 2 L becomes a nonmonotonic function of z. 

We are going to argue that this nonmonotonic behavior is a finite size effect. We 
make the scaling ansatz 

(\n\y\ 2 ) z , L = -a { ^\z)\nL + c(z), (129) 

where £ = 1 if z = 0, 1 and ( = 2 otherwise, while c(z) depends on z but not on L. To 
check the L dependence of (In |^I / | 2 ) 2 L in figure 0, we average (In |^I / | 2 ) 2 L over a narrow 
interval of z's for each L. Figure 0(b) shows the z dependence of a (z) (•) and c(z) 
(■) obtained in this way. In addition, a (z) and c(z) calculated for z — 0, 1 without 
averaging over the narrow interval of z's are shown by open circles and open squares, 
respectively. 

We observe that oc Q (z), if calculated by averaging over a finite range of z's, is almost 
constant and close to a = 2.173. In contrast, a (z = 0, 1) ~ 2.09, if calculated without 
averaging over a finite range of z's, is close to a^ 2 ^ = 2.091. We also find that \c(z)\ 
increases near the boundaries. We conclude that it is the nonmonotonic dependence of 
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\c(z)\ on z that gives rise to the nonmonotonic dependence of (ln|\l/| 2 ) zi on z. This 
finite-size effect is of order 1/ InL and vanishes in the limit L — > oo. 

Boundary condition changing operator 

Next, we impose mixed boundary conditions by either (i) coupling the Z 2 network 
model to an external reservoir through point contacts or (ii) by introducing a long- 
range lead between two nodes from the Z 2 network model, as shown in figure |||. In 
this way, when X > X { , a single Kramers' pair of helical edge states indicated by the 
wavy lines in figure [8] is present on segments of the boundary, while the complementary 
segments of the boundary are devoid of any helical edge state (the straight lines in 
figure H). The helical edge states either escape the Z 2 network model at the nodes at 
which leads to a reservoir are attached [the green lines in figure |8|(a) and figure @(b)], 
or shortcut a segment of the boundary through a nonlocal connection between the two 
nodes located at the corners [figure §](c)]. These are the only options that accommodate 
mixed boundary conditions and are permitted by time-reversal symmetry. As shown 
by Cardy in P5[ , mixed boundary conditions are implemented by boundary-condition- 
changing operators in conformal field theory. Hence, the geometries of figure |^ offer 
yet another venue to test the hypothesis of two-dimensional conformal invariance at the 
metal-insulator quantum critical point. 

When coupling the Z 2 network model to an external reservoir, we shall consider 
two cases shown in figure |8|(a) and figure ||(b), respectively. 

First, we consider the case of figure |](a) in which only two nodes from the Z 2 
network model couple to the reservoir. At each of these point contacts, the scattering 
matrix S that relates incoming to outgoing waves from and to the reservoir is a 2 x 2 
matrix which is invariant under time reversal ^S**^ = S, and it must be proportional 
to the unit 2x2 matrix up to an overall (random) phase. Hence, the two-point-contact 
conductance in the geometry of figure |8](a) is unity, however far the two point contacts 
are from each other. 

Second, we consider the case of figure §(b) in which there are again two point 
contacts, however each lead between the Z 2 network model and the reservoir now 
supports two instead of one Kramers' doublets. The two point-contact scattering 
matrices connecting the Z 2 network model to the reservoirs are now 4x4 matrices, 
which leads to a non-vanishing probability of backscattering. Hence, this even channel 
two-point-contact conductance is expected to decay as a function of the separation 
between the two attachment points of the leads to the network. 

To test whether the two-point-contact conductance in figure |](a) and figure ||(b) 
do differ as dramatically as anticipated, we have computed numerically the two-point- 
contact conductance at the quantum critical point X = X l for a single realization of 
the static disorder. The two-point-contact conductance is calculated by solving for the 



stationary solution of the time-evolution operator with input and output leads [36 



We choose the cylindrical geometry imposed by periodic boundary conditions along the 
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Figure 8. (a) The system with two point contacts (green curves) attached (a) at the 
two interfaces between different types of boundaries and (b) at the reflecting boundary. 
The periodic boundary conditions are imposed for the horizontal direction. The thick 
wavy and solid lines on the edges represent two different types of boundaries, with and 
without a helical edge mode at X > Xi , respectively, (c) Closed network with mixed 
boundaries. Each dashed line or curve represents a Kramers' doublet. 



horizontal directions in figures [|(a) and |^(b) for a squared network with the linear size 
L = 200. Figure ||(a) shows with the symbol • the dependence on r, the distance between 
the two contacts in figure |8](a), of the dimensionless two-point-contact conductance g. 
It is evidently r independent and unity, as expected. Figure |^(a) also shows with 
the symbol o the dependence on r of the dimensionless two-point conductance g for 
leads supporting two Kramers' doublet as depicted in figure §](b). Although it is not 
possible to establish a monotonous decay of the two-point-contact conductance for a 
single realization of the static disorder, its strong fluctuations as r is varied are consistent 
with this claim. 

We turn our attention to the closed geometry shown in figure §|(c). We recall 
that it is expected on general grounds that the moments of the two-point conductance 
in a network model at criticality, when the point contacts are far apart, decay 
as power laws with scaling exponents proportional to the scaling exponents Ag 
f3G| . |37| . Consequently, after tuning the Z 2 network model to criticality, the anomalous 
dimensions at the node (corner) where the boundary condition is changed must vanish, 

Af z 2l°) = 0, (130) 

since the two-point-contact conductance in figure |8](a) is r independent. ( p.30|) is another 
signature of the nontrivial topological nature of the insulating side at the Anderson 
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Figure 9. (a) The distance r dependence of the two-point-contact conductance g for 
the mixed boundary (•) and the reflecting boundary (o). (b) The zero dimensional 
anomalous dimension Aq°' Z2 '°' 1 obtained from the wave function amplitude on the link 
which connecting the boundary condition changing points. 



transition that we want to test numerically. Thus, we consider the geometry figure 
|8](c) and compute numerically the corner anomalous dimensions. This is done using the 
amplitudes of the stationary wave function restricted to the links connecting the two 
corners where the boundary conditions are changed. Figure |9|(b) shows the numerical 
value of the corner anomalous dimension /\^°' Z2 '°\ The linear sizes of the network are 
L = 50, 80, 120, 150, and 180 and the number of disorder realizations is 10 5 for each L. 
We observe that /s!^' 1 " 2 ^ i s zero within the error bars, thereby confirming the validity 
of the prediction (|130Q at the metal-to-Z 2 -topological-insulator transition. 



5. Conclusions 



In summary, we have mapped the Z 2 network model to a 4 x 4 Dirac Hamiltonian. 
In the clean limit of this Dirac Hamiltonian, we expressed the Kane-Mele Z 2 invariant 
as an SU(2) Wilson loop and computed it explicitly. In the presence of weak time- 
reversal symmetric disorder, the NLSM that can be derived out of this Dirac Hamiltonian 
describes the metal-insulator transition in the Z 2 network model and yields bulk scaling 
exponents that belong to the standard two-dimensional symplectic universality class; 
an expectation confirmed by the numerics in [|1^] and |[20|| . A sensitivity to the Z 2 
topological nature of the insulating state can only be found by probing the boundaries, 
which we did numerically in the Z 2 network model by improving the quality of the 
numerical study of the boundary multifractality in the Z 2 network model. 
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